Acanthophora spicifera GROWTH RATE and APICES Analysis, Script Chunks, and Plots

This is the analysis of two runs of Acanthophora spicifera salinity and nutrient experiments conducted on the lanai in St. John 616 in November and December 2023. These experiments incorporated two salinity/nutrient levels and two temperature levels.

Packages loaded:

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(strengejacke)
## # Attaching packages (red = needs update)
## ⚠ ggeffects  1.3.4    ✔ sjlabelled 1.2.0 
## ⚠ sjstats    0.18.2   ✔ esc        0.5.1 
## 
## Warnings or errors in CRAN checks for package(s) 'ggeffects'.
## Update packages in red with 'sj_update()'.
#library(mmtable2) no longer available
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(piecewiseSEM)
library(easystats)
library(magrittr)

Open weight dataset and make columns for growth rate from initial and final weights

acan_growth <- read.csv("../input/acanthophora_growth_r1.csv")

Make a new column for weight change (difference final from initial)

acan_growth$growth_rate_percent <- (acan_growth$final_weight - acan_growth$initial_weight) / acan_growth$initial_weight * 100

Make a new column for daily growth rate from 8 day study (steady growth rate assumed rather than exponential)

acan_growth$steady_growth_daily <- acan_growth$growth_rate_percent / 8

Assigns temperature as a factor

acan_growth$temperature <- as.factor(acan_growth$temperature)

Assigns treatment as characters from integers then to factors

acan_growth$treatment <- as.factor(as.character(acan_growth$treatment))

Assign run as a factor

acan_growth$run <- as.factor(acan_growth$run)

Assign plant ID as a factor

acan_growth$plant_ID <- as.factor(acan_growth$plant_ID)

Assign RLC order as a factor

acan_growth$rlc_order <- as.factor(acan_growth$rlc_order)

Create subsets for the plots and remove one datapoint that is strongly skewing data (overall result does not change)

acan_growth <- subset(acan_growth, secondary_apices != "200")

acan_growth$treatment_graph[acan_growth$treatment == 0] <- "1) 35ppt/0.5umol"
acan_growth$treatment_graph[acan_growth$treatment == 1] <- "2) 28ppt/14umol" 

acan_growth$temperature_graph[acan_growth$temperature == 24] <- "24°C"
acan_growth$temperature_graph[acan_growth$temperature == 28] <- "28°C"

Make a histogram of the growth rate data

hist(acan_growth$growth_rate_percent, main = paste("Acanthophora spicifera 9-Day Growth (%)"), col = "maroon", labels = TRUE)

#or
acan_growth %>% ggplot(aes(growth_rate_percent)) +
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

##Run model without RLC_order as this has little effect, add initial number of apices

growth_model_acan <- lmer(formula = growth_rate_percent ~ treatment + temperature +
                            (1 | plant_ID) + (1 | run) + (1 | initial_axes), data = acan_growth, REML = TRUE)

Plot histograms and performance checks

hist(resid(growth_model_acan))

plot(resid(growth_model_acan) ~ fitted(growth_model_acan))

qqnorm(resid(growth_model_acan))
qqline(resid(growth_model_acan))

performance::check_model(growth_model_acan)

Get R2 and summarize

rsquared(growth_model_acan)
##              Response   family     link method    Marginal Conditional
## 1 growth_rate_percent gaussian identity   none 0.008527329   0.7167619
summary(growth_model_acan)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) +  
##     (1 | run) + (1 | initial_axes)
##    Data: acan_growth
## 
## REML criterion at convergence: 772.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15114 -0.44427  0.05451  0.46396  1.98607 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  plant_ID     (Intercept) 148.58   12.190  
##  initial_axes (Intercept)  49.05    7.003  
##  run          (Intercept)  67.40    8.210  
##  Residual                 105.99   10.295  
## Number of obs: 95, groups:  plant_ID, 48; initial_axes, 8; run, 2
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)
## (Intercept)     -2.864      7.179  1.765  -0.399    0.733
## treatment1       1.229      2.141 44.627   0.574    0.569
## temperature28   -3.348      4.155 38.429  -0.806    0.425
## 
## Correlation of Fixed Effects:
##             (Intr) trtmn1
## treatment1  -0.141       
## temperatr28 -0.288 -0.007

View random effects levels and make a table of the model output

ranef(growth_model_acan)
## $plant_ID
##     (Intercept)
## 1    9.53417840
## 2   11.50686729
## 3    3.87228371
## 4    7.48459427
## 5   -8.96283234
## 6   20.49190957
## 7  -15.56038211
## 8  -12.17546306
## 9   10.63486896
## 10  15.98775355
## 11   6.44751682
## 12  -0.91797790
## 13   4.67935909
## 14   6.86528363
## 15 -11.18300107
## 16   6.05450444
## 17 -10.95411555
## 18  13.84194148
## 19 -18.33904191
## 20  -7.22510672
## 21   8.01159244
## 22 -18.95843963
## 23  -7.19580500
## 24  -1.89213433
## 25  -9.34308528
## 26   0.86169829
## 27  11.39784645
## 28 -11.62101683
## 29   0.73408357
## 30   0.08916723
## 31 -10.19593676
## 32  10.63710508
## 33   6.49379722
## 34  12.18624772
## 35  -0.15971334
## 36  -9.20885772
## 37  10.61902370
## 38   3.56617631
## 39   7.37701722
## 40  -4.37985585
## 41  -9.74270835
## 42  -6.44915512
## 43  -5.04609703
## 44 -21.24271841
## 45  -1.29978106
## 46  -1.83057799
## 47   2.66228410
## 48  11.84670284
## 
## $initial_axes
##    (Intercept)
## 3    3.0315165
## 4    0.7459840
## 5    2.3403631
## 6    0.7300085
## 7    7.3280180
## 8  -11.2173697
## 10   0.2887361
## 11  -3.2472564
## 
## $run
##   (Intercept)
## 1    5.465689
## 2   -5.465689
## 
## with conditional variances for "plant_ID" "initial_axes" "run"
tab_model(growth_model_acan, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  growth_rate_percent
Predictors Estimates std. Error CI Statistic p df
(Intercept) -2.86 7.18 -17.13 – 11.40 -0.40 0.691 88.00
treatment [1] 1.23 2.14 -3.03 – 5.48 0.57 0.567 88.00
temperature [28] -3.35 4.15 -11.60 – 4.91 -0.81 0.423 88.00
Random Effects
σ2 105.99
τ00 plant_ID 148.58
τ00 initial_axes 49.05
τ00 run 67.40
ICC 0.71
N plant_ID 48
N run 2
N initial_axes 8
Observations 95
Marginal R2 / Conditional R2 0.009 / 0.717

Plot the model effects

plot(allEffects(growth_model_acan))

Construct null model to perform likelihood ratio test REML must be FALSE

acan_growth_treatment_null <- lmer(formula = growth_rate_percent ~ temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
acan_growth_model2 <- lmer(formula = growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
anova(acan_growth_treatment_null, acan_growth_model2)
## Data: acan_growth
## Models:
## acan_growth_treatment_null: growth_rate_percent ~ temperature + (1 | plant_ID) + (1 | run)
## acan_growth_model2: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## acan_growth_treatment_null    5 798.17 810.94 -394.08   788.17          
## acan_growth_model2            6 799.41 814.74 -393.71   787.41 0.7559  1
##                            Pr(>Chisq)
## acan_growth_treatment_null           
## acan_growth_model2             0.3846
acan_growth_temperature_null <- lmer(formula = growth_rate_percent ~ treatment + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
acan_growth_model3 <- lmer(formula = growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
anova(acan_growth_temperature_null, acan_growth_model3)
## Data: acan_growth
## Models:
## acan_growth_temperature_null: growth_rate_percent ~ treatment + (1 | plant_ID) + (1 | run)
## acan_growth_model3: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run)
##                              npar    AIC    BIC  logLik deviance  Chisq Df
## acan_growth_temperature_null    5 798.02 810.79 -394.01   788.02          
## acan_growth_model3              6 799.41 814.74 -393.71   787.41 0.6119  1
##                              Pr(>Chisq)
## acan_growth_temperature_null           
## acan_growth_model3               0.4341

Simple boxplots for treatment and temperature

acan_growth %>% ggplot(aes(treatment_graph, growth_rate_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = temperature), position = "jitter", show.legend = TRUE) + 
  labs(x="treatment", y= "9-Day Growth (%)", title= "A", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("T0) 35 ppt/0.5 μmol N", "T1) 28 ppt/14 μmol N")) + 
  ylim(-60, 60) + stat_mean() + 
  scale_color_manual(values = c("azure4", "darkgoldenrod1")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

acan_growth %>% ggplot(aes(temperature_graph, growth_rate_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) + 
  labs(x="Temperature", y= "9-day Growth (%)", title= "B", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("24°C", "28°C")) + 
  ylim(-60, 50) + stat_mean() + 
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  scale_color_manual(values = c("azure3", "darkgoldenrod")) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means

acan_growth %>% group_by(treatment) %>% summarise_at(vars(growth_rate_percent), list(mean = mean))
## # A tibble: 2 × 2
##   treatment  mean
##   <fct>     <dbl>
## 1 0         -2.96
## 2 1         -1.33
acan_growth %>% group_by(temperature) %>% summarise_at(vars(growth_rate_percent), list(mean = mean))
## # A tibble: 2 × 2
##   temperature   mean
##   <fct>        <dbl>
## 1 24          -0.543
## 2 28          -3.73

##RUN MODEL FOR SECONDARY APICES

First plot histogram

acan_growth %>% ggplot(aes(secondary_apices)) +
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Run model for secondary apices

apices_model_acan <- lmer(formula = secondary_apices ~ treatment + temperature +
                            (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = TRUE)

Check performance of model

hist(resid(apices_model_acan))

plot(resid(apices_model_acan) ~ fitted(apices_model_acan))

qqnorm(resid(apices_model_acan))
qqline(resid(apices_model_acan))

performance::check_model(apices_model_acan)

R squared

rsquared(apices_model_acan)
##           Response   family     link method  Marginal Conditional
## 1 secondary_apices gaussian identity   none 0.2748104   0.6747653
summary(apices_model_acan)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: secondary_apices ~ treatment + temperature + (1 | plant_ID) +  
##     (1 | run) + (1 | rlc_order)
##    Data: acan_growth
## 
## REML criterion at convergence: 777
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94937 -0.57432  0.00079  0.45150  2.21087 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  plant_ID  (Intercept)  52.30    7.232  
##  rlc_order (Intercept)  76.24    8.732  
##  run       (Intercept)  49.05    7.003  
##  Residual              144.41   12.017  
## Number of obs: 95, groups:  plant_ID, 48; rlc_order, 24; run, 2
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     40.566      6.385  2.219   6.353 0.018356 *  
## treatment1      21.965      4.338 12.841   5.063 0.000226 ***
## temperature28   -1.007      4.814 22.713  -0.209 0.836247    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmn1
## treatment1  -0.336       
## temperatr28 -0.377 -0.005

View random effects levels, print table and plot all effects

ranef(apices_model_acan)
## $plant_ID
##    (Intercept)
## 1    2.1006896
## 2   -0.4197587
## 3    0.1660575
## 4    0.1660575
## 5    0.4998379
## 6   -3.7009093
## 7    1.5579201
## 8   -3.2729392
## 9   -2.6842007
## 10   3.4068828
## 11   4.1584404
## 12   4.7885525
## 13   4.1142738
## 14  -1.6234315
## 15  -3.4195436
## 16  -0.8990953
## 17  -3.7248190
## 18   2.9963766
## 19   1.2410021
## 20   2.0811516
## 21  -1.3884242
## 22  -6.8493956
## 23  -2.1572453
## 24  -2.1572453
## 25  -5.8447058
## 26  -5.2145937
## 27  -5.4689270
## 28  -5.4689270
## 29   0.7458996
## 30   5.1566842
## 31  -5.9674006
## 32   0.7537950
## 33   2.3927203
## 34   1.7626082
## 35   3.7743900
## 36 -10.2981132
## 37  11.2248719
## 38   0.7230038
## 39   2.2874895
## 40   3.1276389
## 41   5.7628866
## 42   3.2424382
## 43  -6.0742812
## 44   7.1580725
## 45   6.6290198
## 46  -3.2427362
## 47  -3.1714078
## 48   1.0293394
## 
## $rlc_order
##    (Intercept)
## 1   -2.5208706
## 2   -5.7058278
## 3  -11.1497242
## 4   -9.7538717
## 5   -0.8423165
## 6   -2.9409828
## 7    4.7802381
## 8   -7.1586863
## 9    7.9545288
## 10  -0.9602800
## 11  -0.8439845
## 12   4.4926149
## 13   6.6672924
## 14   5.8770567
## 15  14.3796372
## 16  -4.2787335
## 17   2.6274081
## 18  -2.3146177
## 19   9.4375799
## 20   8.7370447
## 21  -2.1057634
## 22  -5.7356272
## 23  -4.9661900
## 24  -3.6759247
## 
## $run
##   (Intercept)
## 1   -4.707122
## 2    4.707122
## 
## with conditional variances for "plant_ID" "rlc_order" "run"
tab_model(apices_model_acan, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  secondary_apices
Predictors Estimates std. Error CI Statistic p df
(Intercept) 40.57 6.39 27.88 – 53.25 6.35 <0.001 88.00
treatment [1] 21.96 4.34 13.34 – 30.59 5.06 <0.001 88.00
temperature [28] -1.01 4.81 -10.57 – 8.56 -0.21 0.835 88.00
Random Effects
σ2 144.41
τ00 plant_ID 52.30
τ00 rlc_order 76.24
τ00 run 49.05
ICC 0.55
N plant_ID 48
N run 2
N rlc_order 24
Observations 95
Marginal R2 / Conditional R2 0.275 / 0.675 
plot(allEffects(apices_model_acan))

Construct null model to perform likelihood ratio test REML must be FALSE

acan_apices_treatment_null <- lmer(formula = secondary_apices~ temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
acan_apices_model2 <- lmer(formula = secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
anova(acan_apices_treatment_null, acan_apices_model2)
## Data: acan_growth
## Models:
## acan_apices_treatment_null: secondary_apices ~ temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## acan_apices_model2: secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## acan_apices_treatment_null    6 821.05 836.37 -404.52   809.05          
## acan_apices_model2            7 805.52 823.40 -395.76   791.52 17.525  1
##                            Pr(>Chisq)    
## acan_apices_treatment_null               
## acan_apices_model2          2.835e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acan_apices_temperature_null <- lmer(formula = secondary_apices ~ treatment + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
acan_apices_model3 <- lmer(formula = secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
anova(acan_apices_temperature_null, acan_apices_model3)
## Data: acan_growth
## Models:
## acan_apices_temperature_null: secondary_apices ~ treatment + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## acan_apices_model3: secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
##                              npar    AIC    BIC  logLik deviance  Chisq Df
## acan_apices_temperature_null    6 803.57 818.89 -395.78   791.57          
## acan_apices_model3              7 805.52 823.40 -395.76   791.52 0.0469  1
##                              Pr(>Chisq)
## acan_apices_temperature_null           
## acan_apices_model3               0.8286

Plots

acan_growth %>% ggplot(aes(treatment_graph, secondary_apices)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = temperature), position = "jitter", show.legend = TRUE) + 
  labs(x="treatment", y= "Number of Secondary Apices", title= "A", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("T0) 35 ppt/0.5 μmol N", "T1) 28 ppt/14 μmol N")) + 
  ylim(-1, 130) + stat_mean() + 
  scale_color_manual(values = c("azure4", "darkgoldenrod1")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

acan_growth %>% ggplot(aes(temperature_graph, secondary_apices)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) + 
  labs(x="Temperature", y= "Number of Secondary Apices", title= "B", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("24°C", "28°C")) + 
  ylim(-1, 130) + stat_mean() + 
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  scale_color_manual(values = c("azure3", "darkgoldenrod")) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means

acan_growth %>% group_by(treatment) %>% summarise_at(vars(secondary_apices), list(mean = mean))
## # A tibble: 2 × 2
##   treatment  mean
##   <fct>     <dbl>
## 1 0          40.1
## 2 1          61.7
acan_growth %>% group_by(temperature) %>% summarise_at(vars(secondary_apices), list(mean = mean))
## # A tibble: 2 × 2
##   temperature  mean
##   <fct>       <dbl>
## 1 24           51.0
## 2 28           50.5